#
# ==============
# Figure A1: Parallel trends Pre-post plot. Relate incidence of
# agency-level hearing to trend of program IP rates.
# ==============
#

  # Program-level data.
  HT = data.table(read_dta(file.path(DATALOC,"ProgramLevel.dta")))
  print(HT)
  HT[,outlaysamountm := as.numeric(outlaysamountm)]
  
  # Number of hearings per agency.
  HT[FY < 2022,num_hearings := sum(agen_nhear_all,na.rm=T),by="agencyid"]
  HT[FY < 2022,num_hearing_years := sum((agen_nhear_all>0),na.rm=T),by="agencyid"]
  # Total outlays across time by program.
  HT[,tot_outlays := sum(outlaysamountm),by="programid"]

  # Year of first hearing per agency.
  HT = merge(HT,HT[(agen_nhear_all>0),.(FY_of_first=min(FY)),by="agencyid"],by="agencyid",all.x=T)
  setkey(HT,agencyid,programid,FY)
  print(HT)

  # Years from first hearing.
  HT[,years_since_first := FY - FY_of_first]
  print(HT)

  # Plots.
  f = file.path(FIGURELOC,"FigureA01.pdf")
  pdf(f,width=10,height=6)
  par.old = par(mar=c(3.1,4.1,0.1,0.1),cex.lab=1.1,cex.axis=1.1)

  # All programs.
  HT[,plot(x=years_since_first,y=iprate,type="n",ann=F,axes=F,log="y")]
  y_at = c(.001,.01,.1,1,10,100)
  axis(1,family="mono");axis(2,at=y_at,labels=y_at,las=2,family="mono")
  abline(v=0,lwd=2,col="gray")
  # Less than 1mm outlays.
  HT[tot_outlays < 1e06,lines(x=years_since_first,y=iprate,type="l",lwd=1,pch=19,cex=.8,col="gray"),by="programid"]
  # More than 1mm outlays.
  HT[tot_outlays > 1e06,lines(x=years_since_first,y=iprate,type="l",lwd=3,pch=19),by="programid"]
  # Unweighted average.
  s = HT[,loess.smooth(x=years_since_first,y=iprate,span=1/3) ]
  lines(x=s$x,y=s$y,lwd=5,col=1)
  title(ylab="Program IP rate (percent, log scale)")
  title(xlab="Years from first hearing",line=2)
  legend("bottomright",lwd=c(1,3,5),col=c("gray","black",1),title="Total outlays",legend=c("Less than $1mm","Greater than $1mm","Unweighted loess"))

  par(par.old)
  dev.off()
